using DataFrames, Plots, Statistics, JSON, Clp
plotlyjs()
include(joinpath(dirname(pwd()),"src/TuLiPa.jl")); # using Dates, JuMP, HiGHS, CSV, Clustering
The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.
Demo 1 gave an introduction to the TuLiPa framework and how a dataset can be built, and how a problem can be built, updated, solved and results queried. Demo 2 will show that TuLiPa also can be used on a more complete dataset, and with more complex model objects.
1.1-1.3 Build, solve and look at solution
1.4 Look at time usage with different types of horizons.
1.5 Compare time usage between JuMP_Prob and HiGHS_Prob.
1.6 Compare time usage with different solvers.
1.7 Generate price-series for hydropower scheduling
This is a simplified version of the dataset NVE uses for its long-term power market analyses (for example https://www.nve.no/energi/analyser-og-statistikk/langsiktig-kraftmarkedsanalyse/). The dataset consist of:
We cannot publish the full dataset, but many of the profiles for wind, solar, hydro and demand can be downloaded from https://www.nve.no/energi/analyser-og-statistikk/vaerdatasett-for-kraftsystemmodellene/.
# The hydropower storages in the dataset needs boundary conditions for the state variables
function addStartEqualStopAllStorages!(modelobjects)
for obj in values(modelobjects)
if obj isa BaseStorage
trait = StartEqualStop(obj)
modelobjects[getid(trait)] = trait
end
end
end
# Power balances needs slack variable for when the inelastic supply (wind, solar, RoR) is higher than the inelastic demand
function addPowerUpperSlack!(modelobjects) # add after object manipulation
for obj in values(modelobjects)
if obj isa BaseBalance
if getid(getcommodity(obj)) == Id("Commodity", "Power")
balancename = getinstancename(getid(obj))
varname = "SlackVar_" * balancename
varkey = Id(FLOW_CONCEPT, varname)
var = BaseFlow(varkey)
sethorizon!(var, gethorizon(obj))
setlb!(var, LowerZeroCapacity())
arrowname = "SlackArrow_" * balancename
arrowkey = Id(ARROW_CONCEPT, arrowname)
arrow = BaseArrow(arrowkey, obj, BaseConversion(PlusOneParam()), 0)
addarrow!(var, arrow)
modelobjects[varkey] = var
end
end
end
end
# Remove start-up costs. Does not make sense to have them when the horizon does not have a fine time-resolution
function remove_startupcosts!(modelobjects)
for (id,obj) in modelobjects
if obj isa StartUpCost
delete!(modelobjects, id)
end
end
end
remove_startupcosts! (generic function with 1 method)
The function runs several scenarios and plot prices, production, demand and hydro storage for each scenario. The function only plots results for chosen price areas.
Unlike Demo 1, the function now support more types of objects. This includes AdaptiveHorizons (horizon built based on residual load), Horizons with offset, SegmentedArrow (for hydropower plants with PQ-curves), BaseAggSupplyCurve (for aggregated power plants) and ExogenBalance (exogen price areas)
# Make list of scenarios
function getscenarios(dt; years)
[TwoTime(getisoyearstart(dt), getisoyearstart(yr)) for yr in years]
end
# Run scenarios and plot results from the whole problem
function runscenarios(scenarios, modelobjects)
runscenarios(scenarios, modelobjects, values(modelobjects))
end
# Run scenarios and plot results from chosen model objects (e.g. only one price area)
function runscenarios(scenarios, modelobjects, resultobjects)
numperiods_powerhorizon = getnumperiods(power_horizon)
numperiods_hydrohorizon = getnumperiods(hydro_horizon)
# We use the datatime for plotting results
dt = getdatatime(scenarios[1])
# Time vector for power_horizon, can be Adaptive
if power_horizon isa AdaptiveHorizon
numperiods_data = Int(getduration(hydro_horizon)/power_horizon.unit_duration)
x1 = [dt + power_horizon.unit_duration*(t-1) for t in 1:numperiods_data]
else
numperiods_data = numperiods_powerhorizon
x1 = [dt + getstartduration(power_horizon, t) for t in 1:numperiods_powerhorizon]
end
# Time vector for hydro_horizon
x2 = [dt + getstartduration(hydro_horizon, t) for t in 1:numperiods_hydrohorizon]
# Horizons can have offset, assumes its the same for both horizons
if hasoffset(hydro_horizon)
x1 = x1 .+ getduration(getoffset(power_horizon))
x2 = x2 .+ getduration(getoffset(hydro_horizon))
end
# Problem can be a HiGHS_Prob or a JuMP_Prob
# mymodel = JuMP.Model(HiGHS.Optimizer)
# set_silent(mymodel)
# prob = JuMP_Prob(collect(values(modelobjects)), mymodel)
@time prob = HiGHS_Prob(collect(values(modelobjects)))
# Order result objects into lists
powerbalances = []
rhsterms = []
plants = []
plantbalances = []
plantarrows = Dict()
demands= []
demandbalances = []
demandarrows = Dict()
hydrostorages = []
for obj in resultobjects # regret using Ids rather than the whole object
# Powerbalances
if obj isa BaseBalance
if getinstancename(getid(getcommodity(obj))) == "Power"
if !isexogen(obj)
push!(powerbalances, getid(obj))
for rhsterm in getrhsterms(obj)
push!(rhsterms,getid(rhsterm))
end
end
end
end
# Hydrostorages
if obj isa BaseStorage
if getinstancename(getid(getcommodity(getbalance(obj)))) == "Hydro"
push!(hydrostorages,getid(obj))
end
end
# Supply and demands
if obj isa BaseFlow
# The type of supply or demand can be found based on the arrows
arrows = getarrows(obj)
# Simple supplies and demands
powerarrowbool = [(getid(getcommodity(getbalance(arrow))) == Id("Commodity", "Power")) & !(arrow isa SegmentedArrow) for arrow in arrows]
powerarrows = arrows[powerarrowbool]
if sum(powerarrowbool) == 1
if isingoing(powerarrows[1])
push!(plants,getid(obj))
push!(plantbalances,getid(getbalance(powerarrows[1])))
elseif !isingoing(powerarrows[1])
push!(demands,getid(obj))
push!(demandbalances,getid(getbalance(powerarrows[1])))
end
end
# Transmissions
if sum(powerarrowbool) == 2
for arrow in arrows
balance = getbalance(arrow)
if getid(getcommodity(balance)) == Id("Commodity", "Power")
if isingoing(arrow) && (balance in resultobjects)
push!(plants,getid(obj))
push!(plantbalances,getid(balance))
if isexogen(balance)
plantarrows[getid(obj)] = arrow
end
elseif !isingoing(arrow) && (balance in resultobjects)
push!(demands,getid(obj))
push!(demandbalances,getid(balance))
if isexogen(balance)
demandarrows[getid(obj)] = arrow
end
end
end
end
end
# Supplies with SegmentedArrows (hydropower with PQ-kurves)
pqarrowbool = [arrow isa SegmentedArrow for arrow in arrows]
pqarrows = arrows[pqarrowbool]
if sum(pqarrowbool) == 1
if isingoing(pqarrows[1])
push!(plants,getid(obj))
push!(plantbalances,getid(getbalance(pqarrows[1])))
end
end
end
# Aggregated supplies (thermal power plants aggregated into one or more equivalent supplies)
# TODO: Result should be a sum of all clusters, not separated
if obj isa BaseAggSupplyCurve
instance = getinstancename(getid(obj))
concept = getconceptname(getid(obj))
balance = getbalance(obj)
for c in 1:getnumclusters(obj)
newname = string(instance,"_",c)
push!(plants,Id(concept,newname))
push!(plantbalances,getid(balance))
end
end
end
# Matrices to store results per time period, scenario and object
prices = zeros(numperiods_data, length(scenarios), length(powerbalances))
rhstermvalues = zeros(numperiods_data, length(scenarios), length(rhsterms))
production = zeros(numperiods_data, length(scenarios), length(plants))
consumption = zeros(numperiods_data, length(scenarios), length(demands))
hydrolevels = zeros(numperiods_hydrohorizon, length(scenarios), length(hydrostorages))
# Update and solve scenarios, and collect results
for (s, t) in enumerate(scenarios)
@time update!(prob, t)
@time solve!(prob)
println("Objective value in scenario $(s): ", getobjectivevalue(prob))
# Collect results for power production/demand/transmission in GW
if power_horizon isa SequentialHorizon
for j in 1:numperiods_powerhorizon
# Timefactor transform results from GWh to GW/h regardless of horizon period durations
timefactor = getduration(gettimedelta(power_horizon, j))/Millisecond(3600000)
# For powerbalances collect prices and rhsterms (like inelastic demand, wind, solar and RoR)
for i in 1:length(powerbalances)
prices[j, s, i] = getcondual(prob, powerbalances[i], j)/1000 # from €/GWh to €/MWh
for k in 1:length(rhsterms)
if hasrhsterm(prob, powerbalances[i], rhsterms[k], j)
rhstermvalues[j, s, k] = getrhsterm(prob, powerbalances[i], rhsterms[k], j)/timefactor
end
end
end
# Collect production of all plants
for i in 1:length(plants) # TODO: Balance and variable can have different horizons
concept = getconceptname(plants[i])
if concept != AGGSUPPLYCURVE_CONCEPT
arrows = getarrows(modelobjects[plants[i]])
pqarrowbool = [arrow isa SegmentedArrow for arrow in arrows]
pqarrows = arrows[pqarrowbool]
if sum(pqarrowbool) == 1
production[j, s, i] = 0
for k in 1:length(getconversions(pqarrows[1]))
segmentid = getsegmentid(pqarrows[1], k)
production[j, s, i] += getvarvalue(prob, segmentid, j)*abs(getconcoeff(prob, plantbalances[i], segmentid, j, j))/timefactor
end
else
if isexogen(modelobjects[plantbalances[i]])
# TODO: Balance and variable can have different horizons
arrow = plantarrows[plants[i]]
horizon = gethorizon(arrow)
conversionparam = _getcontributionparam(arrow)
querytime = getstarttime(horizon, j, t)
querydelta = gettimedelta(horizon, j)
conversionvalue = getparamvalue(conversionparam, querytime, querydelta)
production[j, s, i] = getvarvalue(prob, plants[i], j)*conversionvalue/timefactor
else
production[j, s, i] = getvarvalue(prob, plants[i], j)*abs(getconcoeff(prob, plantbalances[i], plants[i], j, j))/timefactor
end
end
else
production[j, s, i] = getvarvalue(prob, plants[i], j)*abs(getconcoeff(prob, plantbalances[i], plants[i], j, j))/timefactor
end
end
# Collect demand of all demands
for i in 1:length(demands) # TODO: Balance and variable can have different horizons
if isexogen(modelobjects[demandbalances[i]])
arrow = demandarrows[demands[i]]
horizon = gethorizon(arrow)
conversionparam = _getcontributionparam(arrow)
querytime = getstarttime(horizon, j, t)
querydelta = gettimedelta(horizon, j)
conversionvalue = getparamvalue(conversionparam, querytime, querydelta)
consumption[j, s, i] = getvarvalue(prob, demands[i], j)*conversionvalue/timefactor
else
consumption[j, s, i] = getvarvalue(prob, demands[i], j)*abs(getconcoeff(prob, demandbalances[i], demands[i], j, j))/timefactor
end
end
end
else # if power_horizon isa AdaptiveHorizon it is more complex (also test switching order, could be faster to start with objects and traits?)
# Non-consecutive time units have been aggregated in the model, and values need to be distributed to the correct time period
acc = 0
prev_acc = 0
for mp in 1:numperiods_hydrohorizon
for block in 1:power_horizon.num_block
j = (mp-1)*power_horizon.num_block + block
period = power_horizon.periods[j]
timefactor = getduration(period)/Millisecond(3600000)
acc += getlength(period)
for segment in period.units
for unit in segment
fineperiod = prev_acc + unit
# For powerbalances collect prices and rhsterms (like inelastic demand, wind, solar and RoR)
for i in 1:length(powerbalances)
prices[fineperiod, s, i] = getcondual(prob, powerbalances[i], j)/1000
for k in 1:length(rhsterms)
if hasrhsterm(prob, powerbalances[i], rhsterms[k], j)
rhstermvalues[fineperiod, s, k] = getrhsterm(prob, powerbalances[i], rhsterms[k], j)/timefactor
end
end
end
# Collect production of all plants
for i in 1:length(plants) # TODO: Balance and variable can have different horizons
concept = getconceptname(plants[i])
if concept != AGGSUPPLYCURVE_CONCEPT
arrows = getarrows(modelobjects[plants[i]])
pqarrowbool = [arrow isa SegmentedArrow for arrow in arrows]
pqarrows = arrows[pqarrowbool]
if sum(pqarrowbool) == 1
for k in 1:length(getconversions(pqarrows[1]))
segmentid = getsegmentid(pqarrows[1], k)
production[fineperiod, s, i] += getvarvalue(prob, segmentid, j)*abs(getconcoeff(prob, plantbalances[i], segmentid, j, j))/timefactor
end
else
if isexogen(modelobjects[plantbalances[i]])
# TODO: Balance and variable can have different horizons
arrow = plantarrows[plants[i]]
horizon = gethorizon(arrow)
conversionparam = _getcontributionparam(arrow)
querytime = getstarttime(horizon, j, t)
querydelta = gettimedelta(horizon, j)
conversionvalue = getparamvalue(conversionparam, querytime, querydelta)
production[j, s, i] = getvarvalue(prob, plants[i], j)*conversionvalue/timefactor
else
production[fineperiod, s, i] = getvarvalue(prob, plants[i], j)*abs(getconcoeff(prob, plantbalances[i], plants[i], j, j))/timefactor
end
end
else
production[fineperiod, s, i] = getvarvalue(prob, plants[i], j)*abs(getconcoeff(prob, plantbalances[i], plants[i], j, j))/timefactor
end
end
# Collect demand of all demands
for i in 1:length(demands) # TODO: Balance and variable can have different horizons
if isexogen(modelobjects[demandbalances[i]])
arrow = demandarrows[demands[i]]
horizon = gethorizon(arrow)
conversionparam = _getcontributionparam(arrow)
querytime = getstarttime(horizon, j, t)
querydelta = gettimedelta(horizon, j)
conversionvalue = getparamvalue(conversionparam, querytime, querydelta)
consumption[j, s, i] = getvarvalue(prob, demands[i], j)*conversionvalue/timefactor
else
consumption[fineperiod, s, i] = getvarvalue(prob, demands[i], j)*abs(getconcoeff(prob, demandbalances[i], demands[i], j, j))/timefactor
end
end
end
end
end
prev_acc = acc
end
end
# Collect hydro storage levels
for j in 1:numperiods_hydrohorizon
for i in 1:length(hydrostorages)
hydrolevels[j, s, i] = getvarvalue(prob, hydrostorages[i], j)/1000 # Gm3 TODO: convert to TWh with global energy equivalents of each storage
end
end
end
# Only keep rhsterms that have at least one value (TODO: Do the same for sypply and demands)
rhstermtotals = dropdims(sum(rhstermvalues,dims=(1,2)),dims=(1,2))
rhstermsupplyidx = []
rhstermdemandidx = []
for k in 1:length(rhsterms)
if rhstermtotals[k] < 0
push!(rhstermsupplyidx, k)
elseif rhstermtotals[k] > 0
push!(rhstermdemandidx, k)
end
end
# Put rhsterms together with supplies and demands
rhstermsupplyvalues = rhstermvalues[:,:,rhstermsupplyidx]*-1
rhstermdemandvalues = rhstermvalues[:,:,rhstermdemandidx]
rhstermsupplynames = [getinstancename(rhsterm) for rhsterm in rhsterms[rhstermsupplyidx]]
rhstermdemandnames = [getinstancename(rhsterm) for rhsterm in rhsterms[rhstermdemandidx]]
supplynames = [[getinstancename(plant) for plant in plants];rhstermsupplynames]
supplyvalues = reshape([vcat(production...);vcat(rhstermsupplyvalues...)],(numperiods_data,length(scenarios),length(supplynames)))
demandnames = [[getinstancename(demand) for demand in demands];rhstermdemandnames]
demandvalues = reshape([vcat(consumption...);vcat(rhstermdemandvalues...)],(numperiods_data,length(scenarios),length(demandnames)))
# Prepare for plotting results
hydronames = [getinstancename(hydro) for hydro in hydrostorages]
powerbalancenames = [split(getinstancename(powerbalance), "PowerBalance_")[2] for powerbalance in powerbalances]
# Plot results for each scenario
for (s, t) in enumerate(scenarios)
scenyear = string(getisoyear(getscenariotime(t)))
datayear = string(getisoyear(getdatatime(t)))
# Plot prices
display(plot(x1, sum(prices[:,s,:],dims=3), labels=reshape(powerbalancenames,1,length(powerbalancenames)), size=(800,500), title="Prices in " * datayear * " for weather scenario " * scenyear, ylabel="€/MWh"))
# Plot supplies and demands
# supplychart = areaplot(x1, sum(supplyvalues[:,s,:],dims=3),labels=reshape(supplynames,1,length(supplynames)),title="Supply", ylabel = "GWh/h")
# demandchart = areaplot(x1, sum(demandvalues[:,s,:],dims=3),labels=reshape(demandnames,1,length(demandnames)),title="Demand", ylabel = "GWh/h")
# display(plot([supplychart,demandchart]...,layout=(1,2),size=(1600,500)))
supplychart = areaplot(x1, sum(supplyvalues[:,s,:],dims=(3,2)),title="Supply", ylabel = "GWh/h")
demandchart = areaplot(x1, sum(demandvalues[:,s,:],dims=(3,2)),title="Demand", ylabel = "GWh/h")
display(plot([supplychart,demandchart]...,layout=(1,2),size=(900,500)))
# Plot storages (only TWh because of our input data)
display(areaplot(x2, sum(hydrolevels[:,s,:],dims=3),labels=reshape(hydronames,1,length(hydronames)),size=(800,500),title="Reservoir levels", ylabel = "TWh")) #
# Plot list of yearly mean production and demand for each supply/demand
meandemand = dropdims(mean(demandvalues[:,s,:],dims=1),dims=1)
meanproduction = dropdims(mean(supplyvalues[:,s,:],dims=1),dims=1)
supplydf = sort(DataFrame(Supplyname = supplynames, Yearly_supply_TWh = meanproduction*8.76),[:Yearly_supply_TWh], rev = true)
demanddf = sort(DataFrame(Demandname = demandnames, Yearly_demand_TWh = meandemand*8.76),[:Yearly_demand_TWh], rev = true)
supplydf[!,:ID] = collect(1:length(supplynames))
demanddf[!,:ID] = collect(1:length(demandnames))
joineddf = select!(outerjoin(supplydf,demanddf;on=:ID),Not(:ID))
# show(joineddf,allcols=true, allrows=true)
# Check that total supply equals total demand
show(combine(joineddf, [:Yearly_supply_TWh, :Yearly_demand_TWh] .=> sum∘skipmissing))
end
end;
sti_themadata = "data_fra_thema"
tst = JSON.parsefile(joinpath(sti_themadata, "dataset_thema.json"))
themastructure = getelements(tst)
tse = JSON.parsefile(joinpath(sti_themadata, "tidsserier_thema.json"))
themaseries = getelements(tse, sti_themadata)
elements = vcat(themaseries,themastructure)
length(elements)
9426
This would increase the problem size a lot, so we could not continue with this dataset without aggregating it somehow. But it gives the user the possibility to aggregate the hydropower however they want.
NB! We don't continue with this dataset.
sti_themadata = "data_fra_thema"
tst = JSON.parsefile(joinpath(sti_themadata, "dataset_thema_excl_hydro_nose.json")) # European dataset without aggregated hydro in Norway and Sweden
themastructure = getelements(tst)
tse = JSON.parsefile(joinpath(sti_themadata, "tidsserier_thema.json"))
themaseries = getelements(tse, sti_themadata)
dse = JSON.parsefile("data_fra_dynmodell/tidsserier_detd.json")
detdseries = getelements(dse);
dst = JSON.parsefile("data_fra_dynmodell/dataset_detd.json")
detdstructure = getelements(dst);
elements1 = vcat(themaseries,themastructure,detdseries,detdstructure)
length(elements1)
35314
# Select the weather scenarios 1981-2010
# Used by RotatingTimeVector to reduce the amount of data
scenarioyearstart = 1981
scenarioyearstop = 2011
push!(elements, getelement(TIMEPERIOD_CONCEPT, "ScenarioTimePeriod", "ScenarioTimePeriod",
("Start", getisoyearstart(scenarioyearstart)), ("Stop", getisoyearstart(scenarioyearstop))))
# We choose the horizon (time-resolution) of the commodities. We set the duration of the horizons to 1 year
# Having different horizons for different commodities gives flexibility. Here we use it to have more details in the
# variables connected to power, than variables connected to hydro. We also need the storage
# The hydro horizon is sequential with a weekly time-resolution.
hydro_horizon = SequentialHorizon(13*4, Hour(168))
# The power horizon is an AdaptiveHorizon which is dynamically built based on the residual load (see src/horizons.jl for description)
# There are 7 clusters, which mean that the power horizon will have 7 periods inside each hydro horizon period
# For each hydro horizon period, units (hours) are grouped together with units that have a similar residual load, into 7 clusters
rhsdata = DynamicRHSAHData("Power") # alternative StaticRHSAHData("Power", 2021, scenarioyearstart, scenarioyearstop)
# rhsdata = StaticRHSAHData("Power", 2025, scenarioyearstart, scenarioyearstop)
method = KMeansAHMethod() # alternative PercentilesAHMethod()
clusters = 7
unit_duration = Hour(3)
power_horizon = AdaptiveHorizon(clusters, unit_duration, rhsdata, method, 13*4, Hour(168)) #, offset=MsTimeDelta(Week(1)))
# power_horizon = SequentialHorizon(13*4*5, Hour(168))
# Insert horizons into commodities. E.g. all batteries will have the power horizon, since they interact with the power market
function set_horizon!(elements, commodity, horizon)
# If element already exist, replace horizon with new
for element in elements
if element.typename == "BaseCommodity"
if element.instancename == commodity
element.value[HORIZON_CONCEPT] = horizon
return
end
end
end
# Else, add commodity to element list
push!(elements, getelement(COMMODITY_CONCEPT, "BaseCommodity", commodity,
(HORIZON_CONCEPT, horizon)))
end
# The power horizon will be assigned to everything contributing to the power market, or batteries
set_horizon!(elements, "Power", power_horizon)
set_horizon!(elements, "Battery", power_horizon)
# The hydro horizon will be assigned to storages, bypasses and spill variables of hydropower plants (not release because it contributes to the power market and therefore needs to have power_horizon as its horizon)
set_horizon!(elements, "Hydro", hydro_horizon);
modelobjects = getmodelobjects(elements);
length(modelobjects)
1407
When the model objects have been created we can manipulate them however we want. In this example we simplify the problem by aggregating areas and power plants, and removing short term storage systems and start-up costs. We also add boundary conditions to storages.
function simplify!(modelobjects)
# Aggregate price areas and add power balance slack variables
# For the new area FRACHE, the transmission line FRA-CHE is transformed into a demand based on the loss and utilization of the line
aggzoneareadict = Dict("NLDBEL" => ["NLD","HUB_NLD","BEL","HUB_BEL"],
"FRACHE" => ["FRA","CHE"],
"AUTCZE" => ["AUT","CZE"],
"BAL" => ["LTU","LVA","EST","HUB_OST"],
"DMK" => ["DK1","HUB_DK1","DK2","HUB_DK2"],
"NOS" => ["NO1","NO2","NO5"],
"NON" => ["NO3","NO4"],
"SEN" => ["SE1","SE2"],
"SES" => ["SE3","SE4"])
aggzonedict = Dict()
for (k,v) in aggzoneareadict
aggzonedict[Id(BALANCE_CONCEPT,"PowerBalance_" * k)] = [modelobjects[Id(BALANCE_CONCEPT,"PowerBalance_" * vv)] for vv in v]
end
aggzone!(modelobjects, aggzonedict)
# Add slack variable for excessive renewable power
addPowerUpperSlack!(modelobjects)
# Start-up-costs are not compatible with aggregatesupplycurve! or AdaptiveHorizon
remove_startupcosts!(modelobjects)
# Aggregate all simple plants (only connected to power market, mostly thermal) for each area into 4 equivalent plants
aggregatesupplycurve!(modelobjects, 4)
# Short-term storage systems are only needed when the horizon is fine
removestoragesystems!(modelobjects, Hour(10))
# Only calculate AdaptiveHorizon based on residual loads in these areas
residualloadareas!(modelobjects, ["DEU","NLDBEL","POL","GBR","BAL"])
# Storages have state-dependant variables that need a boundary conditions
# We set the starting storage to be equal to the ending storage, x[0] = x[T] (for horizon where t in 1:T)
addStartEqualStopAllStorages!(modelobjects)
end
simplify!(modelobjects);
Due to AdaptiveHorizon we get hourly price volatility even though the LP-problem only has seven periods per week. This leads to seven price levels per week, where hours with similar residual loads will have the same price, because they where grouped together into one period.
# Choose scenarios
scenarios = getscenarios(2025; years=1981:1981)
# Choose areas to see results from
resultobjects = getpowerobjects(modelobjects,["NOS","DMK","DEU"]) #,"NOS","DMK","DEU", "FRACHE"
runscenarios(scenarios, modelobjects, resultobjects)
0.162143 seconds (1.00 M allocations: 81.181 MiB) 0.509293 seconds (4.76 M allocations: 162.221 MiB) 0.940007 seconds Objective value in scenario 1: 5.357321527781434e10
1×2 DataFrame Row │ Yearly_supply_TWh_sum_skipmissing Yearly_demand_TWh_sum_skipmissing │ Float64 Float64 ─────┼────────────────────────────────────────────────────────────────────── 1 │ 862.1 862.1
Still the same amount of periods as the adaptive horizon, but the price has a lot less volatility now.
# Function to change power_horizon in elements list and make list of modelobjects
function get_modelobjects!(elements, new_power_horizon, new_hydro_horizon=nothing)
# Change power_horizon in elements list
set_horizon!(elements, "Power", new_power_horizon)
set_horizon!(elements, "Battery", new_power_horizon)
# Set new power_horizon to the global variable
global power_horizon = new_power_horizon
# If hydro horizon is changed
if !isnothing(new_hydro_horizon)
set_horizon!(elements, "Hydro", new_hydro_horizon)
global hydro_horizon = new_hydro_horizon
end
# Create modelobjects and manipulate
modelobjects = getmodelobjects(elements)
simplify!(modelobjects)
return modelobjects
end
power_horizon = SequentialHorizon(52*7, Day(1))
modelobjects = get_modelobjects!(elements, power_horizon)
# Choose areas to see results from
resultobjects = getpowerobjects(modelobjects,["DEU","DMK","NOS"]) #,"NOS","DMK","DEU", "FRACHE"
runscenarios(scenarios, modelobjects, resultobjects);
0.170732 seconds (1.18 M allocations: 86.933 MiB) 0.301693 seconds (3.26 M allocations: 114.350 MiB) 0.941034 seconds Objective value in scenario 1: 5.320472699786423e10
1×2 DataFrame Row │ Yearly_supply_TWh_sum_skipmissing Yearly_demand_TWh_sum_skipmissing │ Float64 Float64 ─────┼────────────────────────────────────────────────────────────────────── 1 │ 863.472 863.472
# Function to change power_horizon in elements list and make list of modelobjects
function get_detailed_modelobjects!(elements, new_power_horizon, new_hydro_horizon)
# Change power_horizon in elements list
set_horizon!(elements, "Hydro", new_hydro_horizon)
set_horizon!(elements, "Power", new_power_horizon)
set_horizon!(elements, "Battery", new_power_horizon)
# Set new power_horizon to the global variable
global hydro_horizon = new_hydro_horizon
global power_horizon = new_power_horizon
# Create modelobjects and manipulate
modelobjects = getmodelobjects(elements)
# Add slack variable for excessive renewable power
addPowerUpperSlack!(modelobjects)
# Storages have state-dependant variables that need a boundary conditions
# We set the starting storage to be equal to the ending storage, x[0] = x[T] (for horizon where t in 1:T)
addStartEqualStopAllStorages!(modelobjects)
return modelobjects
end
hydro_horizon = SequentialHorizon(14, Day(1))
power_horizon = SequentialHorizon(14*24, Hour(1))
modelobjects = get_detailed_modelobjects!(elements, power_horizon, hydro_horizon)
# Choose areas to see results from
resultobjects = getpowerobjects(modelobjects,["DEU","DK1","NO2"]) #,"NOS","DMK","DEU", "FRACHE"
runscenarios(scenarios, modelobjects, resultobjects);
2.137076 seconds (4.66 M allocations: 392.566 MiB, 52.44% gc time) 0.189255 seconds (2.53 M allocations: 72.127 MiB) 16.106352 seconds Objective value in scenario 1: 2.3719842714987435e9
1×2 DataFrame Row │ Yearly_supply_TWh_sum_skipmissing Yearly_demand_TWh_sum_skipmissing │ Float64 Float64 ─────┼────────────────────────────────────────────────────────────────────── 1 │ 860.172 860.172
Sequential daily is fastest since it updates the problem with time-series data for daily periods.
Sequential hourly is slowest since it updates the problem with time-series data for hourly periods. It would give the best prices, but at the the cost of being too slow.
Adaptive takes a bit more time to update than sequential daily since we read and cluster hourly time-series data into seven periods per week. However the solvetime is similar to Sequential daily. This shows that Adaptive could be a good compromise between quality and speed.
function time_usage_horizon(modelobject, scenarios)
@time prob = HiGHS_Prob(collect(values(modelobjects)))
for (s, t) in enumerate(scenarios)
@time update!(prob, t)
@time solve!(prob)
end
end
new_hydro_horizon = SequentialHorizon(52, Hour(168))
display("Sequential daily")
new_power_horizon = SequentialHorizon(52*7, Day(1))
modelobjects = get_modelobjects!(elements, new_power_horizon, new_hydro_horizon)
time_usage_horizon(modelobjects, scenarios)
display("Sequential hourly")
new_power_horizon = SequentialHorizon(52*7*24, Hour(1))
modelobjects = get_modelobjects!(elements, new_power_horizon)
time_usage_horizon(modelobjects, scenarios)
display("Adaptive cluster hours into 7 periods per week")
new_power_horizon = AdaptiveHorizon(clusters, unit_duration, DynamicRHSAHData("Power"), KMeansAHMethod(), 52, Hour(168))
modelobjects = get_modelobjects!(elements, new_power_horizon)
time_usage_horizon(modelobjects, scenarios)
"Sequential daily"
0.174841 seconds (1.18 M allocations: 86.933 MiB) 0.274568 seconds (3.26 M allocations: 114.350 MiB) 0.954871 seconds
"Sequential hourly"
6.430420 seconds (35.37 M allocations: 2.084 GiB, 20.40% gc time) 10.504482 seconds (90.18 M allocations: 2.807 GiB, 31.21% gc time) 69.449921 seconds
"Adaptive cluster hours into 7 periods per week"
0.151033 seconds (1.00 M allocations: 81.181 MiB) 0.512200 seconds (4.76 M allocations: 162.221 MiB) 0.932611 seconds
We made HiGHS_Prob because HiGHS was the fastest open-source solver possible to use with JuMP, but was very slow to update the problem in the JuMP implementation. We therefore implemented HiGHS_Prob directly on the HiGHS package. (Some new PRs has partly solved this problem, but HiGHS_Prob is still noticeably faster.)
We time building, updating and solving the problem for different Probs and solvers, to see the difference. This is of course only a small problem, and solving times could be improved by choosing other solver settings.
function time_usage_probs(modelobjects, scenarios, opt)
if opt == "H"
display("HiGHS_Prob")
@time prob = HiGHS_Prob(collect(values(modelobjects)))
elseif opt == "JC"
display("JuMP_Prob w Clp")
model = Model(Clp.Optimizer)
set_silent(model)
@time prob = JuMP_Prob(collect(values(modelobjects)), model)
elseif opt == "JH"
display("JuMP_Prob w HiGHS")
model = Model(HiGHS.Optimizer)
set_silent(model)
@time prob = JuMP_Prob(collect(values(modelobjects)), model)
end
for (s, t) in enumerate(scenarios)
@time update!(prob, t)
@time solve!(prob)
end
end
time_usage_probs(modelobjects, scenarios, "H")
time_usage_probs(modelobjects, scenarios, "JH")
time_usage_probs(modelobjects, scenarios, "JC")
"HiGHS_Prob"
0.162058 seconds (1.00 M allocations: 81.181 MiB) 0.505826 seconds (4.76 M allocations: 162.221 MiB) 0.939175 seconds
"JuMP_Prob w HiGHS"
0.289726 seconds (3.11 M allocations: 128.728 MiB) 0.970802 seconds (5.12 M allocations: 166.085 MiB) 1.566489 seconds (2.32 M allocations: 99.554 MiB)
"JuMP_Prob w Clp"
0.277916 seconds (2.80 M allocations: 122.394 MiB) 0.978563 seconds (5.12 M allocations: 166.085 MiB) 2.607031 seconds (467.22 k allocations: 102.040 MiB)
The results show that HiGHS_Prob is fastest for building, updating and solving the problem. This is because it is custom-built for our purpose. The higher time usage of JuMP_Prob w HiGHS actually comes from transferring the JuMP model into HiGHS rather that solving the problem with HiGHS.
We also see that HiGHS_Prob takes more time updating the problem than solving it. This could be improved by using the time-series more efficiently (now querying from time-series is coded very generic). Some possible improvements:
We compare time usage of HiGHS_Prob and JuMP_Prob with Clp, for different solver implementations: dual simplex and interior-point method (IPM). The result shows that HiGHS_Prob with dual simplex is fastest for this specific problem given the default solver settings. With other setting or a different problem (longer horizon / higher resolution / more objects), other combinations could be the fastest.
Overview of solver settings for Clp and HiGHS.
function time_usage_solver(modelobjects, scenarios, opt)
if opt == "HS"
display("HiGHS_Prob dual simplex")
@time prob = HiGHS_Prob(collect(values(modelobjects)))
Highs_setStringOptionValue(prob, "solver", "simplex")
elseif opt == "HI"
display("HiGHS_Prob ipm")
@time prob = HiGHS_Prob(collect(values(modelobjects)))
Highs_setStringOptionValue(prob, "solver", "ipm")
elseif opt == "JCS"
display("JuMP_Prob w Clp dual simplex")
model = Model(Clp.Optimizer)
set_optimizer_attribute(model, "SolveType", 0)
set_silent(model)
@time prob = JuMP_Prob(collect(values(modelobjects)), model)
elseif opt == "JCI"
display("JuMP_Prob w Clp ipm")
model = Model(Clp.Optimizer)
set_optimizer_attribute(model, "SolveType", 4)
set_silent(model)
@time prob = JuMP_Prob(collect(values(modelobjects)), model)
end
for (s, t) in enumerate(scenarios)
@time update!(prob, t)
@time solve!(prob)
end
end
time_usage_solver(modelobjects, scenarios, "HS")
time_usage_solver(modelobjects, scenarios, "HI")
time_usage_solver(modelobjects, scenarios, "JCS")
time_usage_solver(modelobjects, scenarios, "JCI")
"HiGHS_Prob dual simplex"
0.169661 seconds (1.00 M allocations: 81.181 MiB) 1.258970 seconds (4.76 M allocations: 162.221 MiB, 55.89% gc time) 0.937266 seconds
"HiGHS_Prob ipm"
0.153379 seconds (1.00 M allocations: 81.181 MiB) 0.504131 seconds (4.76 M allocations: 162.221 MiB) 3.800384 seconds
"JuMP_Prob w Clp dual simplex"
0.278241 seconds (2.80 M allocations: 122.394 MiB) 0.992424 seconds (5.12 M allocations: 166.085 MiB) 2.607747 seconds (467.22 k allocations: 102.040 MiB)
"JuMP_Prob w Clp ipm"
0.278965 seconds (2.80 M allocations: 122.394 MiB) 0.995247 seconds (5.12 M allocations: 166.085 MiB) 2.367710 seconds (467.22 k allocations: 102.040 MiB)
In demo 4 we need an exogen price to optimize a watercourse against. We generate a 15 year price series for Denmark using the model year 2025 and the weather scenarios 1981-1995.
function make_price_series(elements)
# Horizons can also be made from a starting and ending DateTime, and the duration D of the periods.
# This will generate a horizon of N periods of duration D, and the last period will be the remaining duration (which is smaller than D)
weeks = 4
hydro_horizon = SequentialHorizon(getisoyearstart(1981), getisoyearstart(1996), Week(weeks))
power_horizon = AdaptiveHorizon(7, Hour(1), DynamicRHSAHData("Power"), KMeansAHMethod(), getisoyearstart(1981), getisoyearstart(1996), Week(weeks))
set_horizon!(elements, "Power", power_horizon)
set_horizon!(elements, "Battery", power_horizon)
set_horizon!(elements, "Hydro", hydro_horizon)
modelobjects = getmodelobjects(elements)
simplify!(modelobjects)
@time prob = HiGHS_Prob(collect(values(modelobjects)))
datayear = getisoyearstart(2025)
scenarioyear = getisoyearstart(1981)
t = FixedDataTwoTime(datayear, scenarioyear)
@time update!(prob, t)
@time solve!(prob)
# We query the price from Denmark (since NOS price is too flat)
id = Id(BALANCE_CONCEPT, "PowerBalance_DMK")
balance = modelobjects[id]
horizon = gethorizon(balance)
# Generate and store price param as elements
_values = [getcondual(prob, id, s) for s in 1:getnumperiods(horizon)]
generate_param(_values, scenarioyear, horizon, id)
end
# Function to generate and store price param as elements
function generate_param(values::Vector{Float64}, dt::DateTime, horizon::AdaptiveHorizon, id::Id)
# Distribute values in AdaptiveHorizon periods to hours that represents the clusters
points = []
for t in 1:getnumperiods(horizon)
value = values[t]
timedelta = gettimedelta(horizon, t)
startduration = getstartduration(horizon, t)
for unit_range in timedelta.units
for i in unit_range
point_dt = dt + startduration + ((i-1) * timedelta.unit_duration)
point = (point_dt, value)
push!(points, point)
end
end
end
sort!(points)
_index = [Dates.format(p[1], "yyyy-mm-dd HH:MM:SS") for p in points] # TODO: Change to Dates default format?
_values = [p[2] for p in points]
# Plot
display(plot(_index, _values, label=getinstancename(id), size=(900,500)))
# Store in dataelements
e = []
push!(e, (TIMEINDEX_CONCEPT,"VectorTimeIndex","IndexPriceDMK",
("Vector",_index)))
push!(e, (TIMEVALUES_CONCEPT,"VectorTimeValues","ValuesPriceDMK",
("Vector",_values)))
push!(e, (TIMEVECTOR_CONCEPT,"RotatingTimeVector","ProfilePriceDMK",
(TIMEINDEX_CONCEPT,"IndexPriceDMK"),(TIMEVALUES_CONCEPT,"ValuesPriceDMK")))
push!(e, (PARAM_CONCEPT,"MeanSeriesParam","PriceDMK",
("Level",1),
("Profile","ProfilePriceDMK")))
open("priceDMK.json","w") do f
JSON.print(f, e)
end
end
make_price_series(elements)
0.688229 seconds (4.33 M allocations: 319.431 MiB) 11.526492 seconds (78.36 M allocations: 1.873 GiB, 8.52% gc time) 6.059711 seconds